library(vegan)
library(tidyverse)
library(reshape2)
library(scales)
library(forcats)
otu <-
    read.table(
        file = "./metaphlan.tsv",
        sep = "\t",
        header = T,
        row.names = 1,
        quote = ""
    )
group <-
    read.table(
        file = "./group_info.tsv",
        sep = "\t",
        header = T,
        # row.names = 1,
        quote = ""
    ) %>% filter(., !grepl("control", Group)) # remove control
# import clinical (without sex and age)
clinic_original <-
    readxl::read_xlsx("./clincal.xlsx",
        sheet = 1
    )
# import sex and age
clinic_sexage <-
    readxl::read_xlsx("./clincal.xlsx",
        sheet = 2
    ) %>%
    filter(!grepl("C.*", m0)) %>%
    pivot_longer(
        cols = 1:3,
        names_to = "Group",
        values_to = "samples",
        values_drop_na = T
    )
path <-
    read.table(
        file = "./path.csv",
        sep = ",",
        header = T,
        row.names = 1,
        quote = ""
    ) %>% t()
# to a list
clinic <- list(
    TG = clinic_original[, 1:2] %>% drop_na(2) %>% mutate(TG_g = ifelse(TG < 1.7, 1, 2) %>%  as.character()),
    Chol = clinic_original[, c(1, 3)] %>% drop_na(2) %>% mutate(Chol_g = ifelse(Chol < 5.2, 1, 2) %>% as.character()),
    HDL = clinic_original[, c(1, 4)] %>% drop_na(2) %>% mutate(HDL_g = ifelse(HDL < 1, 1, 2) %>% as.character()),
    LDL = clinic_original[, c(1, 5)] %>% drop_na(2) %>% mutate(LDL_g = ifelse(LDL < 3.4, 1, 2) %>% as.character()),
    ALB = clinic_original[, c(1, 6)] %>% drop_na(2) %>% mutate(ALB_g = ifelse(ALB < 35, 1, 2) %>% as.character()),
    Cr = clinic_original[, c(1, 7)] %>% drop_na(2) %>% mutate(Cr_g = ifelse((Cr > 106 & sample %>% str_detect(".*M.*")) | (Cr > 97 & sample %>% str_detect(".*F.*")), 2, 1) %>%  as.character()),
    CTnT = clinic_original[, c(1, 8)] %>% drop_na(2) %>% mutate(CTnT_g = ifelse(cTnT < 0.5, 1, 2) %>%  as.character()),
    BNP = clinic_original[, c(1, 9)] %>% drop_na(2) %>% mutate(BNP_g = ifelse(BNP < 100, 1, ifelse(BNP < 400, 2, 3)) %>%  as.character()),
    HbA1 = clinic_original[, c(1, 10)] %>% drop_na(2) %>% mutate(HbA1_g = ifelse(HbA1c < 5.6, 1, 2) %>%  as.character()),
    sex = clinic_sexage[, c(4, 1)] %>% mutate(gender = ifelse(sex == "male", 1, 2) %>%  as.character()) %>% select(c(-2)),
    age = clinic_sexage[, c(4, 2)],
    subject = clinic_original[, 1] %>% mutate(subject = str_sub(sample, 2, 3) %>% as.numeric() %>%  as.character())
)
# select species from otu
tax_s <-
    slice(otu, str_which(rownames(otu), ".*s__.*")) %>%
    select(., contains(group$Sample)) %>%
    t()
# select genus from otu
tax_g <-
    slice(otu, str_which(rownames(otu), ".*g__((?!s__).)*$")) %>%
    select(., contains(group$Sample)) %>%
    t()
# prepare a matrix
adonis_result <-
    matrix(
        nrow = 12,
        ncol = 6,
        dimnames = list(
            names(clinic),
            c(
                "tax_s_R2",
                "tax_s_p",
                "tax_g_R2",
                "tax_g_p",
                "path_R2",
                "path_p"
            )
        )
    )

{
    adonis_result[1, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["TG"]]$sample, ] ~ TG_g,
            clinic[["TG"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[2, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["Chol"]]$sample, ] ~ Chol_g,
            clinic[["Chol"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[3, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["HDL"]]$sample, ] ~ HDL_g,
            clinic[["HDL"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[4, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["LDL"]]$sample, ] ~ LDL_g,
            clinic[["LDL"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[5, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["ALB"]]$sample, ] ~ ALB_g,
            clinic[["ALB"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[6, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["Cr"]]$sample, ] ~ Cr_g,
            clinic[["Cr"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[7, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["CTnT"]]$sample, ] ~ CTnT_g,
            clinic[["CTnT"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[8, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["BNP"]]$sample, ] ~ BNP_g,
            clinic[["BNP"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()

    adonis_result[9, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["HbA1"]]$sample, ] ~ HbA1_g,
            clinic[["HbA1"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[10, 1:2] <-
        adonis2(tax_s ~ gender,
            clinic[["sex"]][clinic[["sex"]]$samples %in% rownames(tax_s), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[11, 1:2] <-
        adonis2(tax_s ~ age,
            clinic[["age"]][clinic[["sex"]]$samples %in% rownames(tax_s), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[12, 1:2] <-
        adonis2(tax_s[rownames(tax_s) %in% clinic[["subject"]]$sample, ] ~ subject,
            clinic[["subject"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
}
{
    adonis_result[1, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["TG"]]$sample, ] ~ TG_g,
            clinic[["TG"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[2, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["Chol"]]$sample, ] ~ Chol_g,
            clinic[["Chol"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[3, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["HDL"]]$sample, ] ~ HDL_g,
            clinic[["HDL"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[4, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["LDL"]]$sample, ] ~ LDL_g,
            clinic[["LDL"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[5, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["ALB"]]$sample, ] ~ ALB_g,
            clinic[["ALB"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[6, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["Cr"]]$sample, ] ~ Cr_g,
            clinic[["Cr"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[7, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["CTnT"]]$sample, ] ~ CTnT_g,
            clinic[["CTnT"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[8, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["BNP"]]$sample, ] ~ BNP_g,
            clinic[["BNP"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()

    adonis_result[9, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["HbA1"]]$sample, ] ~ HbA1_g,
            clinic[["HbA1"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[10, 3:4] <-
        adonis2(tax_g ~ gender,
            clinic[["sex"]][clinic[["sex"]]$samples %in% rownames(tax_g), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[11, 3:4] <-
        adonis2(tax_g ~ age,
            clinic[["age"]][clinic[["sex"]]$samples %in% rownames(tax_g), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[12, 3:4] <-
        adonis2(tax_g[rownames(tax_g) %in% clinic[["subject"]]$sample, ] ~ subject,
            clinic[["subject"]],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
}
{
    adonis_result[1, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["TG"]]$sample, ] ~ TG_g,
            clinic[["TG"]][clinic[["TG"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[2, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["Chol"]]$sample, ] ~ Chol_g,
            clinic[["Chol"]][clinic[["Chol"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[3, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["HDL"]]$sample, ] ~ HDL_g,
            clinic[["HDL"]][clinic[["HDL"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[4, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["LDL"]]$sample, ] ~ LDL_g,
            clinic[["LDL"]][clinic[["LDL"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[5, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["ALB"]]$sample, ] ~ ALB_g,
            clinic[["ALB"]][clinic[["ALB"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[6, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["Cr"]]$sample, ] ~ Cr_g,
            clinic[["Cr"]][clinic[["Cr"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[7, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["CTnT"]]$sample, ] ~ CTnT_g,
            clinic[["CTnT"]][clinic[["CTnT"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[8, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["BNP"]]$sample, ] ~ BNP_g,
            clinic[["BNP"]][clinic[["BNP"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()

    adonis_result[9, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["HbA1"]]$sample, ] ~ HbA1_g,
            clinic[["HbA1"]][clinic[["HbA1"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[10, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["sex"]]$samples, ] ~ gender,
            clinic[["sex"]][clinic[["sex"]]$samples %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[11, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["age"]]$samples, ] ~ age,
            clinic[["age"]][clinic[["sex"]]$samples %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
    adonis_result[12, 5:6] <-
        adonis2(path[rownames(path) %in% clinic[["subject"]]$sample, ] ~ subject,
            clinic[["subject"]][clinic[["subject"]]$sample %in% rownames(path), ],
            permutations = 999,
            distance = "bray"
        )[1, c(3, 5)] %>% as.matrix()
}
# R2 value
adonis_R2 <- melt(adonis_result[, c(1, 3, 5)])
# create a fake R2 value matrix to adjust color
adonis_R2_f <- melt(adonis_result[, c(1, 3, 5)]) %>%  mutate(Var1=factor(Var1))
adonis_R2_f$value[adonis_R2$value > 0.04] <-
    adonis_R2$value[adonis_R2$value > 0.04] %>% rescale(to = c(0.04, 0.05)) # 让颜色分布均匀
# p value
adonis_p <- melt(adonis_result[, c(2, 4, 6)])
adonis_p$value <-
    ifelse(adonis_p$value <= 0.05,
        ifelse(
            adonis_p$value <= 0.01,
            ifelse(adonis_p$value <= 0.001, "***", "**"),
            "*"
        ),
        NA
    )
# join R2 and p value tables
adonis_p$Var2 <- str_replace_all(adonis_p$Var2, "_p", "_R2")
adonis_R2 <-
    left_join(adonis_R2,
        adonis_p,
        by = c("Var1", "Var2"),
        suffix = c("_R", "_p")
    ) %>%  mutate(Var1=factor(Var1))

p <- ggplot(adonis_R2_f, aes(x = Var2, y = Var1 , fill = value)) +
    geom_tile() +
    geom_text(
        aes(label = label_percent(accuracy = 0.1)(round(
            adonis_R2$value_R, 3
        ))),
        color =
            ifelse(adonis_R2$value_R > 0.025, "white", "black"),
        size = 4
    ) +
    geom_text(
        aes(label = adonis_R2$value_p),
        color =
            ifelse(adonis_R2$value_R > 0.025, "white", "black"),
        size = 4.5,
        vjust = -0.2
    ) +
    scale_fill_distiller(direction = 1) +
    scale_y_discrete(position = "right",limits=factor(adonis_R2$Var1)) +
    coord_fixed() +
    theme_minimal() +
    theme(legend.position = "none")+
    labs(x=NULL, y=NULL)
ggsave(p,filename= "perm.pdf")
# adon <- function(x) {
#   temp <- matrix(ncol = 2,dimnames = list(NA,c("R2","p")))
#   # temp_name <- c(TG,Chol,HDL,LDL,ALB,Cr,CTnT,BNP,EF,HbA1c)
#   for (i in 1:length(clinic)) {
# print(i)
#     temp[i,] <-
#       adonis2(x ~ c(TG,Chol,HDL,LDL,ALB,Cr,CTnT,BNP,EF,HbA1c)[i],
#               clinic[[i]],
#               permutations = 999,
#               distance = 'bray')
#   }
#   temp
# }
# sapply(c(tax_s,tax_g), adon)

